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Abstract. The influence of nucleons superfluidity on the beta relaxation time of degenerate neutron star cores, 
composed of neutrons, protons and electrons, is investigated. We numerically calculate the implied reduction 
factors for both direct and modified Urea reactions, with isotropic pairing of protons or anisotropic pairing of 
neutrons. We find that due to the non-zero value of the temperature and/or to the vanishing of anisotropic gaps 
in some directions of the phase-space, superfluidity does not always completely inhibit beta relaxation, allowing 
for some reactions if the superfluid gap amplitude is not too large in respect to both the typical thermal energy 
and the chemical potential mismatch. We even observe that if the ratio between the critical temperature and the 
actual temperature is very small, a suprathermal regime is reached for which superfluidity is almost irrelevant. 
On the contrary, if the gap is large enough, the composition of the nuclear matter can stay frozen for very long 
durations, unless the departure from beta equilibrium is at least as important as the gap amplitude. These results 
are crucial for precise estimation of the superfluidity effect on the cooling/slowing-down of pulsars and we provide 
online subroutines to be implemented in codes for simulating such evolutions. 
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1. Introduction 

In the simplest model, the so-called npe matter model, 
neutron star (NS) cores are transparent for neutrinos 
and mainly consist of degenerate neutrons, with a small 
admixture of equal numbers of protons and electrons. In 
such conditions, the composition is characterized by the 
fraction of protons among the total number of nucleons, 
x p . Due to beta reactions, the value of x p depends on 
density and is usually determined by the condition of 
beta equilibrium, for instance in the calculations of 
non-dynamical quantities, such as the equation of state 
(EOS). However, the assumption of beta equilibrium is 
only valid if the density of a matter element changes on 
a timescale t p much longer than the beta equilibration 
timescale 773. But in strongly degenerate npe matter, Tp is 
made much longer than in normal matter by the decrease 
of the available phase-space due to the Pauli exclusion 
principle. As a consequence, various astrophysical sce- 
narios in which the condition Tp <C t p is violated were 
pointed out by several authors. Such a situation occurs 
in gravitational collapse of neutron star (Haensel [1992 
Gourgoulhon & Haensel 1993), during pulsar spin-down 
(Reisenegger 119951 119971 Fernandez & Reisenegger 2005), 
or in the neutron star interiors just due to the existence of 
even relatively slow hydrodynamic flows (for an example, 
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see Urpin & Shalybkov I1996f) or of millisecond oscil- 
lations (Reisenegger & Goldreich lrM?! Villain et al. l2Htj5|l . 

Nonetheless, beta equilibrium breaking and change 
of composition are crucial to take into account not only 
from those mainly "mechanical" points of view, but 
also in the thermal evolution of not too old neither too 
young NSs. Indeed, the cooling of such NSs proceeds 
through the emission of neutrinos by the cold (i.e., 
degenerate) npe matter, whose deviation from beta 
equilibrium is characterized by the chemical potential 
mismatch 5u = u n — u p — u p . While in beta equilibrium 
only the non-zero value of the temperature T determines 
the size of the available phase-space, and consequently 
all reactions rates, Sp ^ opens additional volume in 
the phase-space for beta processes. This results in an 
increase of the neutrino emissivity, and as beta reactions 
occurring off-equilibrium produce entropy and therefore 
heat neutron star matter, the net effect is the matter 
heating (Haensel 119921 Gourgoulhon & Haensel 1993 
Reisenegger 119951 119971 Fernandez & Reisenegger 2005). 

Furthermore, the description of any basic physical 
process that takes place inside NSs can be made more 
difficult by another phenomenon, which is the possible 
superfluidity of nucleons. Indeed, it is expected that in 
about one year after neutron star birth, neutrons (pro- 
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tons) in the core become superfluid, when temperature 
goes below critical one, T cn (T cp ). Then, energy gaps, 
A„ and A p , appear in the nucleon excitation energy 
spectra, reducing again the available phase-space. This 
strongly slows down the beta processes in dense matter 
for T <C T c and also makes the dynamical properties, 
such as pulsations (Tee 119951 Andersson & Comer 170011 
Prix & Rieutord[2002), more complicated to evaluate. 

To summarize, the beta processes rate, in its full 
complexity, involves at least four parameters with the 
dimension of energy: S/j,, fceT, A„, and A p . The roles of 
the first and the second pairs of parameters are opposite: 
S/i, kftT increase the available phase-space, while the 
superfluid energy gaps inhibit the reactions. In the limit 
of small deviations from the thermodynamic equilibrium, 
Sfi/ksT <C 1, the beta relaxation rate can be linearized 
in this small parameter and the dissipation can be 
represented by the bulk viscosity (Haensel et al. [2000 
2001] [20021 and references therein). Nevertheless, in the 
general case, both the heating rate and neutrino losses 
have to be calculated explicitly in the out of beta equi- 
librium npe matter. For a non-superfluid matter this was 
done in Haensel (1992), Gourgoulhon & Haensel (JT993), 
Reisenegger (|1995p and Fernandez & Reisenegger (J2005 ). 
The effect of nucleon superfluidity on beta relaxation 
rates was studied using a very crude model by Reisenegger 
(1997). He considered deviation from beta equilibrium 
implied by the the pulsar spin-down, and assumed that 
for 5/j, < A p + A n the beta reactions were completely 
blocked (no reactions), while for 5fi > A p + A n the effect 
of superfluidity could be neglected (normal matter rates). 
As we should discuss in more details later, such a step-like 
modeling of the reduction factor behaviour is unrealistic. 

In the present article, we calculate numerically the 
superfluid reduction factors for a broad range of the 
relevant parameters making also available on-line subrou- 
tines to proceed to any farther calculation. Both direct 
Urea and modified Urea processes are considered in the 
presence of various types of superfluidity. In this way, 
we get correct reduction factors that can be used for 
simulating the evolution of superfluid neutron star cores 
which are off beta equilibrium. 

The formulae for the rates of non-equilibrium beta pro- 
cesses in a non-superfluid npe matter are reminded in Sect. 
[3 Basic features of nucleon superfluidity in NSs are pre- 
sented in Sect. [21 In Sect. 01 the formulae for the super- 
fluid reduction factors are derived, with some illustrative 
numerical results in Sect. [SJ In Sect. we discuss our 
results and their possible application in numerical simu- 
lations of neutron star evolution and dynamics. Finally, 
in the Appendix, the subroutines to calculate the reduc- 
tion factors for such numerical simulations are briefly de- 
scribed, together with some comments on their practical 
use and the address of the website from which they can 
be downloaded. 



2. Out-of-equilibrium beta processes without 
superfluidity 

Non-equilibrium beta processes in normal (non super- 
fluid) neutron star cores have already been the subject of 
several studies fHaensel 119921 Reisenegger I199f)l) . In this 
section, we shall not discuss the astrophysical conditions 
in which the breaking of chemical equilibrium can occur, 
but focus on its microphysical description. Moreover, 
we just give here a brief summary of known results, 
introducing our notations in units h = k& = c = 1 (with 
exception of Sect.OJ, mainly following the conventions of 
Haensel (TT3321) . 

We take the basic npe model of nuclear NS matter and 
assume that 

- the core's content is in a stationary state (no oscil- 
lation), with all particles comoving. The outcome of a 
possible difference between the motions (e.g., rotation) 
of n and p fluids (a crucial phenomenon for superfluids) 
will be considered elsewhere; 

- the nucleons are strongly interacting non-relativistic 
Fermi liquids and the electrons form an ultrarelativis- 
tic ideal Fermi gas. The temperature is sufficiently 
low for the npe components to be strongly degenerate, 
and we shall then work for all of them with quanti- 
ties at zero temperature and at thermodynamic, but 
not necessarily chemical, equilibrium: beta equilibrium 
is not assumed. Hence, the Fermi momenta p-pi, for 
the i species, versus the respective densities Uj, are 
p Fi = (37r 2 ni) 1/3 ; 

- local electromagnetic equilibrium is reached and mat- 
ter is electrically neutral; 

- the neutrinos freely escape from the star. 

Cooling of NSs involves various reactions, among which 
the Urea processes are the only ones that we should deal 
with in this article. Indeed, as neutrinos freely escape from 
the star, other reactions, such as bremsstrahlung, do not 
change the chemical composition of npe matter and are 
then not affected by beta equilibrium breaking. Yet, a 
distinction has to be made between two types of Urea 
processes, since the direct Urea reactions (called there- 
after Durca, following K. Levenfish, see, Yakovlev et al. 
2001) are kinematically allowed only if the proton frac- 
tion is high enough. If not, the modified Urea reactions 
(Murca) 1 play the key-role, nonetheless they can be ne- 
glected when Durca occur. From a practical point of view, 
it means that there is a threshold density below which one 
has to deal with (and only with) Murca reactions, while, 
for densities above that value, it is sufficient to consider 
only Durca processes. As the latter are the simplest, we 
shall start with their description. 



with, in respect to Durca, an additional spectator nucleon. 
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2.1. Direct Urea processes 



with the dimcnsionless integral 



In npe matter, Durca reactions are allowed if 
PFn < PFp + PFe- Assuming local neutrality of matter, 
it is equivalent to p-p n < 2 PFp and then to n p > n/9, 
where n is the total baryonic density n = n p + n n . This 
leads to a threshold value of the density that depends 
on the equation of state, and mainly on the symmetry 
energy, but which is typically several times the standard 
nuclear matter density no = 2.8 x 10 14 gcm~ 3 . 

When beta equilibrium is assumed, the two Durca re- 
actions, 

D n : n — > p + e~ + v e , 

D p : p + e~ — » n + v e , 



(1) 



have the same rates and form, all together, the usual 
Durca process. As we do not assume beta equilibrium, 
we have to consider two different rates (numbers of reac- 
tions per cm 3 and during one second) (see Haensel [1992 
Reisenegger 1 1 9 9 5 ) 



= 



dp e dp p dp n dp v 



(2tt) 3 (2tt) 3 (2tt) 3 (2tt) 3 
x (27r)H(E f -E i )S(P f 



(l-/e)(l -/„)/„ 
Pr)\M lf \ 2 (2) 



and 



Fd„ — 



dp e dp p dp n dp v 



(1 fn)fpfe 



(2tt) 3 (2tt) 3 (2tt) 3 (2tt) 3 
x (2n)U(E f -E l )S(P f -P t )\M lf \ 2 . (3) 

In these equations, |-Mj/| 2 is the squared transi- 
tion amplitude summed and averaged on spins states, 
fi = (l + e^ ei ~^ i ^ T )~ 1 is the Fermi-Dirac distribution 
for a fermion i with energy Si, chemical potential /Xj, and 
temperature T, while the 5-functions ensure the conser- 
vation of the total energy and momentum. The physical 
quantity we will focus on is ATp, that we define as the 
difference between T u n and T d p ■ Ar d=Td„ — ^d p ■ 
This quantity characterizes the relaxation rate (and hence 
timescale) for reaching beta equilibrium. 

In the context of strongly degenerate matter, there is 
a standard approximation, that can be applied here, to 
go beyond the previous formulae in an analytical way, 
the so-called phase-space decomposition (see Shapiro & 
Teukolskv 1 1 983fl . It mainly consists in replacing in the in- 
tegrals all smooth functions by their values at pi = p-pi, 
which enables to factorise out from the integrals the mean 
value of the microphysical factor | A/^ / 1 2 - We shall not 
make explicit the whole calculation, and send the reader 
to Yakovlev et al. H2001[l for a detailed review. However, 
with the following dimensionless variables (i = n,p, e) 



Mi 



where we defined S{i = fi n — ftp ~ Me, we finally get 
Ar D (£) - AT Da I D (0, 



(4) 



OO +OG 



{fn (1 - fp) (1 - fe) $ (X n - Xp - X e - X u + C) 

- fp fe (1 ~ fn) $ {X n + X v - X p - X e + £)} . (6) 

The rate ArD , a physical factor that depends on 
the EOS, takes such typical values that the relax- 
ation time, once the integral © calculated, is around 
- 20 T 9 " 4 s, with T g = T/(10 9 K) (Yakovlev et al. 

EHH). 

As far as the calculation of the integral is con- 
cerned, using {(xi) = fi = 1 — {(—Xi) combined with 
known results on Fermi integrals, we are lead to 



Id(0 = J dx v xi {J D (a^ - - Jd(3!« +0} » ( 7 ) 
o 

where Jd(^) = i( x ) (x 2 + 71-2 The latter integration can 
now be done analytically and gives (Reisenegger[l995) 



Id(0 = £^ (1 
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2.2. Modified Urea 



Durca processes are by many orders of magnitude the 
fastest beta processes in NS cores. However, as already 
mentioned, they can be kinematically forbidden. If this is 
the case, the Murca reactions prevail, which involve an 
additional spectator nucleon to ensure the conservation of 
both energy and momentum. This additional nucleon can 
be either a proton or a neutron. Hence, the distinction 
has to be made between two branches, and not only when 
some nucleons are superfluid. The Murca reactions are 
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(9) 



where N is the spectator nucleon. If this particle is 
a proton, these reactions are called proton branch of 
Murca, and they are called neutron branch whether it is 
a neutron 2 . 

If none of the nucleons is superfluid, we do not really 
have to make distinction between the two branches in the 
phase-space integral, and using the same notations as in 
the Durca case (skipping the details of the calculations 
that are very similar), we can write 



Ar M (0 = AiWoMO 



(10) 



2 Notice that the proton branch has a density threshold, but 
it is much smaller than the threshold for Durca, with the con- 
(5) dition n p > n/65 in npe matter (Yakovlev & Levenfish|1995 ). 
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with 



oo +00 

Im{<0 = J dx v xl J dx n dx p dx e dx Nl dx Nf 

-co 

{/n/jVi(l-/ P )(l-/e)(l-/jV/) 
X S (x n + X Nl —Xp — X e — X N f — x u + £) 

~ fp fe fm (1 - fn) (1 - /n/) 

x 5 (x n + x Nf + x v - Xp - x Nl - x e + £)} . (11) 

Here, x Ni arid XNf are respectively the u x variables" , 
defined as in Eq.J3}, for the initial and final spectator 
nucleons. Finally, Ar Mo is a constant that depends on the 
EOS and takes such values that the relaxation timescale 
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Equation ^TJ can also be written 



OO 

(0 = J dx u 



{J M (a:*-0 - Jm(^ + 0) (12) 



where Jm{x) = f(x) (x A 
to fReisenegger I1995jl 



10tt 2 x 2 + 9tt 4 )/24, which leads 



Im(0 



367 g 7T 6 
1512 



1 
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21 £ 4 



3£ 6 



367 7T 2 367 7T 4 1835 tt 6 

(13) 

In the presence of superfluidity, the rates can be much 
more complicated to calculate. Indeed, in all the previous 
calculations, we could quite easily get rid of all angular 
integrals as they always gave contributions that could be 
factorised out. The factorisation was made possible since 
the only function that depends on the respective directions 
of momenta was, a priori, the square of the transition am- 
plitudes. Indeed, it can be shown, cither in the Durca or 
in the Murca case, that it is a proper approximation to 
use some suitably defined angle averaged values (Friman 
& Maxwell rjlZni Yakovlev & Levenfish ITMB) of \M lf \ 2 , 
instead of doing the full angular integral. However, this 
approximation can be no longer sufficient to factorise the 
angular integrals if some nucleons are superfluid. This re- 
sults from the fact that the Cooper pairing can occur in 
anisotropic states, a feature that is explained in the next 
section dealing with nucleon superfluidity in NS cores. 

3. Superfluidity of nucleons 

Superfluidity of nucleons in NS cores is reviewed by 
Lombardo & Schulze (20(nj and we shall only give here 
a brief summary of the results useful in the following. 
Neutrons are believed to form Cooper pairs due to their 
interaction in the triplet 3 P2 state, while protons form 
singlet 1 So pairs. In the study of the triplet-state neutron 
pairing, one should distinguish the cases of the different 
possible projections mj of nn-pair angular momentum 
J onto a quantization axis z (see, e.g., Amundsen and 
0stgaard 1985): \mj\ = 0, 1, 2. The actual (energetically 



Table 1. Three types of superfluidity. From Haensel et 
(2000) with the kind permission of the authors. 



Type 


Pairing state 




fc B T c /A(0) 


A 
B 

C 


3 P 2 (raj = 0) 
3 P 2 (\mj\ = 2) 


1 

(1 + 3 cos 2 tf) 
sin 2 i? 


0.5669 
0.8416 
0.4926 



most favorable) state of nn-pairs is not known, being 
extremely sensitive to the (still unknown) details of nn 
interaction. One cannot exclude that this state varies 
with density and is a superposition of states with different 
mj. 

Hence, in the following, we shall conform ourselves to 
what is usually done in the community (see, e.g., Yakovlev 
et al. I2001f) and deal with three different superfluidity 
types: 1 So, 3 P2 (mj = 0) and 3 P2 (\mj\ — 2), denoted as 
A, B and C, respectively (Tabic [TJ. The superfluidity of 
type A is attributed to protons, while types B and C may 
be attributed to neutrons and are sufficient to have a gen- 
eral idea about possible impact of neutrons superfluidity 
since in these two cases its effect arc qualitatively different. 

In addition, it can be shown (see, e.g., Yakovlev et al. 
2001} that in our context, it is sufficient, from the mi- 
croscopic point of view, to consider as the only effect of 
superfluidity the introduction of an energy gap 6 in mo- 
mentum dependence of the nucleon energy, e(p). Near the 
Fermi level (\p — p-p\ <§C pp), this dependence can be writ- 
ten as 



e = n - ^J5 2 + v f (p-pf) 2 at p <pf , 
e = fx + \J& 2 + Vp(p -pf) 2 at p>pf, 



(14) 



where pf and vf are the Fermi momentum and Fermi 
velocity of the nucleon, respectively, while \i is the nucleon 
chemical potential. One has <5 2 = A 2 (T)F($), where 
A(T) is the part of the gap's amplitude that depends 
on the temperature, and F(i!)) specifies dependence of 
the gap on the angle •& between the particle momentum 
and the z axis (Table QJ. In case A the gap is isotropic, 
and 5 — A(T). In cases B and C, the gap depends on 
Note that in case C the gap vanishes at the poles of 
the Fermi sphere at any temperature: Fc(0) = Fc{it) = 0. 

Notice also that the gap amplitude, A(T), is derived 
from the standard equation of the BCS theory (see, 
e.g., Yakovlev et al. [T999), with the value of A(0) that 
determines the critical temperature T c . The values of 
fcBT c /A(0) for cases A, B and C are given in Tabled 

For further analysis it is convenient to introduce the 
dimcnsionless quantities: 



ACQ 

knT ' 



T 



V 



knT 



(15) 
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The dimensionless gap y can be presented in the form: 
Va = Va, Vb — v B vl + 3 cos 2 ??, ?/c=vcsin?? (16) 

where the dimensionless gap amplitude v depends only 
on r. In case A the quantity v coincides with the isotropic 
dimensionless gap, while in cases B and C it represents, re- 
spectively, the minimum and maximum gap (as a function 
of ■&) on the nucleon Fermi surface. As we shall see later, 
this implies that for global integrated quantities (such as 
rates), with a given value of v, one can expect a stronger 
effect of superfluidity for the case B and a weaker for the 
case C, while the case A should be in between. Notice fi- 
nally that the dependence of v on r was fitted by Levenfish 
& Yakovlcv (1994). However, their fits will not be used in 
the following. Indeed, we shall only deal with microscopic 
calculations and not with astrophysical applications, and 
we shall therefore use v as a free parameter. 

4. Out of equilibrium beta processes with 
superfluidity 

4.1. General features of the effect of superfluidity on 
reactions rates 

The influence of superfluidity on the beta relaxation time 
was already investigated by Reisenegger (|1997J> . but only 
in a very rough way. The principle that he used is based 
on what was explained in Section [21 i.e., on the fact that 
to estimate this influence on averaged physical quantities, 
such as emissivity or reaction rates, it is sufficient, at the 
first level of approximation, to only consider the creation, 
near the Fermi surface, of a gap in the dispersion relation 
£j = £i(p)- This property can be demonstrated using 
Bogoliubov transformations (see, e.g., Yakovlev et al. 
2001) and has, as a first obvious consequence, strong de- 
crease of the available phase-space. Hence, as Reisenegger 
concluded, this can result in a blocking of the reactions 
during quite long times, even if chemical equilibrium is 
broken. Nevertheless, there are two shortcomings in the 
assumption of a completely frozen composition, as made 
by Reisenegger. First, the temperature is not exactly 
null, which implies the existence of excited states and, 
therefore, of some reactions. Second, the Cooper pairing 
does not always occur in a singlet spin state, which has for 
consequence a possible anisotropy of the gap and of the 
Fermi surface, allowing for reactions in some directions of 
the momentum space. 

Thus, to improve the study of Reisenegger (1997J) , we 
shall apply in the following the procedure that consists 
in adding a gap close to the Fermi surface in the disper- 
sion relation, but do the calculations of relaxation rates 
through integrals similar to those described in Sect. [21 In 
practice, it means, for superfluid nucleons, to replace in 
the integrals the usual dimensionless variables Xi by Zi, 
with 

Zi= E% T M * = sign (a;,) yj x\ + yf , (17) 



where 

x ^ v Fi(Pi - PFi ) and ^M^)^ (T)FiW . 

(18) 

In these expressions, 

- VFi = (dsi/ dpi) Pi=PFi is the Fermi velocity, making 
the definitions for Xj, Eq.QJ and Ea. H18l) . equivalent 
near the Fermi surface when there is no gap; 

- 5i(T, $) is the superfluid gap with T the temperature 
and $ the angle between the particle momentum and 
the quantization axis (see Sect. 03); 

- Vi(T) is the dimensionless gap amplitude function and 
Fi(i9) the dimensionless anisotropy function with, for 
the isotropic case (A), F;(i?) = 1 (see Sect.[3J). 

The next issue that needs to be clarified before doing 
any calculation is the identification of the relevant cases 
to deal with. For instance, as it is well-known that 
protons can only pair in spin-singlets, there is no need 
to evaluate the effect of two simultaneous anisotropic 
pairings. Furthermore, a first approximation, which 
should not be too bad, is to assess that there is always 
only one type of superfluid nucleons, which is equivalent 
to assume that the dominant gap prevails. Hence, we 
shall in the following only study cases with superfluidity 
of one type of nucleons, beginning with Durca processes 
that involves less nucleons than Murca processes. Finally, 
notice that even if our goal is to calculate the superfluid 
equivalent of the rates J3J) and (|10|) . we are able, with 
the phase-space approximation, to extract from these 
rates all microphysical uncertainties (ATd and ATm ), 
leaving dimensionless integrals [Id(0 and Im(Q]- As a 
consequence, in order to also minimize, in the superfluid 
cases, the uncertainties coming from microphysics, it is 
more relevant not to directly calculate the superfluid 
rates, but the ratios between their values with and 
without superfluidity. In such a way, only uncertainties 
linked with the gap still affect the phase-space integrals 
and the results. 

This leads us to write, for all type of reactions, equa- 
tions like 

Ar^ = AT Xo f x (Z,vj) = AT X Bk(Z, Vj ) , (19) 

where Rx(^,Vj) = I x (^,Vj)/ lx{0 are the quantities we 
shall estimate: the reduction factors for the X type of 
reactions. In this expression, i labels the type of superflu- 
idity (A, B or C), and j the type of nucleon that is su- 
perfluid (n or p) with a gap Vj [see definition in Ea. (|18|) ]. 
Finally, notice that by definition, for T greater than both 
critical temperatures, one has R x (£, Vj ) = R x (£, 0) = 1, 
and for T smaller than at least one critical temperature 
R^vj) < 1. 

4.2. Durca 

Taking into account the fact that in the Durca case with 
only one type of superfluidity there might be only one 
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type of anisotropy induced by pairing, we can write the 
reduction factor Rd(£,Uj) a s 



1 r dn 



Id(OJ 4vr 

ir/2 



(20) 



1 



dtf sin(0) H* D (6«i)> 



where, for type A of superfluidity, HJ-,(^, Vj) does not 
depend on the t? angle between the momentum of the 
Cooper pair and the quantization axis. 

Following Eq.lJHJl, one then has 



Hb(£j v j) = J x l J Ax n dx p dx e 



x {/„ (1 - fp) (1 - f e ) 6 (z n - z p - x e - x u + £) 

- fp fe (1 - fn) 5 (z n + X v - Z p - X e + £)} , (21) 

where the z variable for the non-superfluid nucleon has to 
be replaced with the corresponding x variable, while Vj 
is "included" in the z variable for the superfluid nucleon 
and f p = f(z p ) = (l + e 2p ) _1 with a similar definition 
for /„. 

The integration over the electron variable can easily 
be done due to the Dirac functions, and then, using the 
well-known formula for Fermi integrals, 



da; i(x) f(y — x) = 



ev - 1 



= G(y), (22) 



enables us to integrate over the x variable for the non- 
superfluid nucleon, giving 



+oo 

dx v x 2 v / dxj 





x {f( Zj ) G{x v - Zj) + i(-Zj) G(x v - £ + Zj) 
- t(zj) G{x v +£- zj) - {(-zj) G(x v + Zj)} . (23) 

This formula is the simplest that can be reached an- 
alytically and will be evaluated with numerical methods 
(see Section EJ. 

4.3. Murca 

Whatever the situation with superfluidity, the main 
change between Durca and Murca reactions is the pres- 
ence of spectator nucleons that make the calculations 
more difficult. Nevertheless, a great difference is that, 
in the superfluid case, the additional dimensions of the 
integral no longer lead to calculations that can quite 
easily be done analytically, mainly due to anisotropies. 



Anyway, with the same method as previously de- 
scribed, Ea. lfm is now replaced with 

AF Mx = AT Mxo P Mx (£, Vj) = AT Mx R^ x (£, Vj ) , (24) 

with R|j x (£,«,•) = lM x ((i,Vj)/lM(€), while X is the 
index of the branch (n or p), i the index of the super- 
fluidity type and j the index of the superfluid nucleon. 
Here again, for T greater than the critical temperature, 
one has R : M (£, Vj) — W M (^,0) = 1 and in other cases 

Rki x (^) < i- 

Notice that P M Vj) has to include some angular 
factor Am xi similar to df2/(47r) in Eg. (|20fl . but somehow 
more complicated. Hence, we have 




dX;y X 



5 +°° 

11 / dxjfj 

J = 1 -oo 



3=1 



where we have introduced a condensed notation for all 
particles, with exception of the neutrinos : j is the index 
of the particles with j = 1 .. 5, respectively, corresponding 

5 

to n,p, Ni, Nf,e. With this notation, Y[ dflj is a global 

i=i 

angular element of integration, while all other notations 
(zj, fj) agree with previously defined notations with 
the convention that for a non-superfluid nucleon or an 
electron Zj = Xj. Finally, notice that the At: factor in 
front of the expression comes from the integration over 
possible directions of the neutrino momentum. 

To complete this section, we will give more advanced 
analytical formulae corresponding to some specific cases 
we will deal with in the following, not only because they 
are the easiest to treat numerically, but also because it is 
sufficient due to the still huge error-bars concerning dense 
matter microphysics. Those cases are 

- the neutron branch with isotropic superfluidity of pro- 
tons (gap A); 

- the proton branch with any superfluidity of neutrons; 

- the neutron branch with isotropic superfluidity of neu- 
trons (gap A). 

In the first case, due to full isotropy, one can write 

+ OG 

I P M„&Vp) = JL j dx p {{z p ) (H(* p + 0-H(z p -0) , 



(26) 



with 



R(x)= dss 2 ; \ - Us~x) 2 +4tt 2 ) (27) 

J exp(s — x) — 1 v ' 

o 



(25) 
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and 



3£ 6 



1835 tt £ 



(28) 



3lU 5 252 ^ 367 7T 2 367 tt 4 
In the second case, one has 

1 +oo 

^(^») = 37T / dc J dx n i(z n ) 



x (H(z„+£)-H(z„-£)) , (29) 

where c= cos(z?) appears in the y n variable while g(£) and 
H(x) are the same as in Eq.(J2EJl. 

The third case leads to an integral to evaluate numeri- 
cally that is more complicated, since analytical integration 
can be performed only on the electron and non-superfluid 
nucleon variables. With the usual technique, one gets 



rnA / £ 
1 M n Is : v n 



1 



+oo 



7m (0 

— oo 

(H(z nl + z n2 + z n3 + £) - H(z nl + z n2 + z n3 - £)) 



with 



H(x) = / ds s 



s — x 



exp(s — x) — 1 



(31) 




dx n i dx n2 dx n3 f(z„i) f(z n2 ) f(z n3 ) 

(30) 



Fig. 1. Three dimensional view of the reduction factor for 
Durca reaction with isotropic superfluidity of the protons. 
Remember that here by T we mean fcb T due to the chosen 
units. 



and where the normalisation function 7m (0 was defined 
in Ea.ipjl. 

For the numerical calculation, as was done for 
Durca reactions [see Ea. (|23|) ]. all those integrals for 
Murca reactions are cut in pieces in such a way that 
the Xi variables take only positive values. However, 
we shall not give here more details on the algo- 
rithms that we used, and we send the reader to 
the Appendix or to the subroutine comments on 
' http:/ /luth2. obspm.fr/~etu/villain/Micro/Reduction. html 
Instead, the next section will now describe the main nu- 
merical results. 



5. Numerical results 

Due to the complexity of the multidimensional integra- 
tions, we decided to test several numerical integration al- 
gorithms (see the Appendix for more details) and finally 
we obtained results in perfect agreement in the physi- 
cal parameter range of interest. Here, we shall just give 
a graphical overview of those results for 

- Durca reaction with superfluidity of each type; 

- Murca reaction with isotropic superfluidity of the non- 

spectator nucleon (e.g. proton for the neutron branch); 

- Murca reaction with isotropic superfluidity of the spec- 

tator nucleon (e.g. neutron for the neutron branch). 



The first graph (Fig. ^) depicts in a three-dimensional 
way the reduction factor for the simplest situation, which 
is Durca reaction with isotropic superfluidity (of protons) . 
We shall not give more 3D figures like this one, because 
they are not the best to look at the influence of the type 
of superfluidity, even if they help to have a quick idea 
of the effect of superfluidity. Indeed, this figure shows 
that for huge departures from equilibrium, superfluidity 
really has an impact only when the amplitude of the gap 
is of the order of Sfi, or is much larger than Sfi. In the 
first case, there is an exponential decay of the reduction 
factor, which leads, for the second case, to the expected 
/'frozen composition" (reduction factor equal to 0). Yet, 
for very small values of the gap's amplitude (in respect to 
the temperature), the ratio with Sfi does not matter and 
superfluidity is not relevant since the reduction factor is 
very close to 1. 

Fig. El shows exactly the same results, but in another 
way. Here, the variable is Sfi, and each curve corresponds 
to a given value of the ratio gap /temperature. For huge 
gaps, Logio(gapfT) > 1, it can be noticed that as soon as 
8 ii > gap, and until the reduction factor becomes close to 
1, the curves are parallel. The meaning of this feature of 
the graphs is that, at least for that part of the curves, some 
functions f and g exist, such that the reduction factor can 
be written as R ~ f(gap) Q(gap/5n) . Furthermore, com- 
parison of Figs. and 0] with Fig.|2]enables to see that the 
effect of anisotropy of the gap cannot easily be evaluated 
since it can either increase (type C) or decrease (type B) 
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Reduction Factor for several values of Log 1() (gap/T) 
Durca Isotropic 




Reduction Factor for several values of Log ]() (gap/T) 
Durca Anisotropic Gap B 
' ' ' ' - 0.5- 




-1 1 



Fig. 2. Reduction factor as a function of the departure 
from beta equilibrium for several values of the gap's am- 
plitude. Durca reaction with isotropic superfluidity of the 
protons. Remember that here by T we mean kh T due to 
the chosen units. 



the reduction factor in respect to its value in the isotropic 
case. However, what seems to be the most important is 
the maximal value of the gap on the Fermi surface since 
the value for case B is smaller than for case A, while for 
case C it is larger. As far as Murca reactions are concerned 
(Figs. |S] and EJ), the situation is the following. If only one 
nucleon (the non-spectator one) is superfluid, then the dif- 
ference with the Durca case seems very small and almost 
not visible (see Figs. [2and[SJ|. However, if the spectator 
nucleons (the most represented) arc also superfluid the im- 
pact is much larger, as is depicted by Fig.^l Nevertheless, 
it should not be forgotten that the reduction factors are 
calculated for relaxation times that are really different : 
Durca reactions allow for a much faster relaxation towards 
equilibrium (see, e.g., sect. 3.5 of Yakovlev et al. 1200 II and 
references therein). 

6. Discussion and conclusions 

We calculated the effect of nucleon superfluidity on the 
rates of beta processes in neutron star matter which is 
off beta equilibrium. Superfluidity reduces these rates by 
a factor that depends on the types of the beta process 
(direct Urea or modified Urea) and of superfluidity, but 
also on the dimensionless gaps and chemical-potential 
mismatch parameters A/ksT and 8^,/k^T, respectively. 
Due to the degeneracy of nuclear matter, the reduction 
factors do not depend explicitly on the EOS (see Sect. 14. 

In order to reduce the number of parameters and 
the technical difficulties in calculating the integrals to 
evaluate, we assumed that we could always consider 
only one type of nucleon being superfluid, the larger gap 
prevailing. Such an approximation, which seems to be 



Fig. 3. Reduction factor as a function of the departure 
from beta equilibrium for several values of the gap's am- 
plitude. Durca reaction with anisotropic superfluidity type 
B. Notation as in Fig.|21 



Reduction Factor for several values of Log ]() (gap/T) 
Durca Anisotropic, Gap C 




Log |0 (5n/T) 



Fig. 4. Reduction factor as a function of the departure 
from beta equilibrium for several values of the gap's am- 
plitude. Durca reaction with anisotropic superfluidity type 
C. Notation as in Fig. [5] 



quite reasonable for the nuclear matter of NS cores (see 
Lombardo & Schulze 2001), indeed allows to analytically 
integrate more variables (see Sect. 0J, strongly reducing 
the computational time needed. 

In the case of the direct Urea process, reduction 
factors values were calculated for the isotropic spin- 
singlet pairing of protons and for two types of the 
angle-dependent gaps resulting from the spin-triplet 3 P2 
pairing of neutrons. Since, in the modified Urea case, 
the phase-space integrals, which can be of up to the 
12 th order, are very complex, we decided to first neglect 
the superfluid gaps anisotropy. This approximation can 
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Reduction Factor for several values of Log 1() (gap/T) 
Murca proton branch, isotropic superfluidity of neutrons 




Fig. 5. Reduction factor as a function of the departure 
from beta equilibrium for several values of the gap's am- 
plitude. Murca reaction with isotropic superfluidity of one 
single nucleon. Notation as in Fig. |2] 



Reduction Factor for several values of Log 1() (gap/T) 
Murca neutron branch, isotropic superfluidity of neutrons 




Fig. 6. Reduction factor as a function of the departure 
from beta equilibrium for several values of the gap's ampli- 
tude. Murca reaction with isotropic superfluidity of three 
nucleons out of four which are involved. Notation as in 
Fig.0 



seems quite rough, however, it is supported by the fact 
that anisotropy of the gap can either increase or decrease 
the effect of superfluidity in respect to the case of the 
isotropic gap (see Sect. EJ, while the actual state of the 
triplet neutron Cooper pairs is still poorly known (see 
Sect. 0. 

Subroutines for the numerical calculation of the 
superfluid reduction factors are briefly described in the 
Appendix with link to the website on which they can be 
found. They yield the reduction factors for given A/fceT 
and 5n/kftT. One has then to multiply by these factors 



the analytical expressions for the non-superfluid beta re- 
action rates to get the physical rates for superfluid matter. 

Our results, combined with an equation of state of the 
neutron star core, can be used in numerical simulations 
of neutron star pulsations, to account for the dissipation 
due to the non-equilibrium beta processes implied by the 
local baryon density variations. Such a modeling will be 
valid both in a highly nonlinear suprathermal regime of 
dissipation Sfi > k^T as well as in the 5fi <C k^T limit 
where the bulk viscosity description can be used (see 
Haensel et al. HUM HUUTl I2UU2|) . Another application 
concerns the calculation of non-equilibrium heating and 
neutrino emission of superfluid matter in a spinning down 
pulsar, which would allow one to model in a realistic 
manner the cooling/slowing-down of superfluid pulsars. 
Both problems will be studied by us in the near future. 

The formulae presented in the present article were ob- 
tained under the simplifying assumption that neutrons, 
protons, and electrons move together as one single fluid, 
i.e., they are all at rest in one single reference system 
(local rest-frame). On a dynamic timescale, such an ap- 
proximation is correct if all three fluids are normal (see, 
e.g., Villain et al. E005). However, in the most gen- 
eral case of superfluid neutron and protons, we have to 
deal with three fluids, and therefore with three differ- 
ent flows of matter. First, we have a normal fluid com- 
posed of electrons and of normal components of neutrons 
and protons (this is a fluid of excited neutron and pro- 
ton quasi-particles). Then, we have two flows of superfluid 
condensates, the neutron and proton one, connected be- 
tween themselves by the superfluid entrainmcnt (Andreev- 
Bashkin effect). If we connect the local rest-frame with 
normal fluid, then the equilibrium distributions of neu- 
tron and proton quasiparticles in the condensates should 
be calculated taking into account the relative velocities 
u ne — u n — u e and u pe = u p — u e , respectively. For ex- 
ample, for isotropic gaps, in the rest-frame comoving with 
normal fluid, the neutron superfluid quasiparticle energies 
aree^(p) = e n (p-mu ne )+^ nn p-u ne +j np p-u pe , where m 
is neutron mass, and terms quadratic in u ne and u pe were 
neglected; the last two terms describe the entrainment ef- 
fects and 7 Q /3 (a, (3 — n,p) is the microscopic entrain- 
mcnt matrix calculated by Gusakov & Haensel (2005). 
Generally, macroscopic flows velocities are small compared 
to the Fermi velocities and lowest-order approximations 
are valid. The calculations in the three-fluid model will be 
the next step of our study of non-equilibrium beta pro- 
cesses in superfluid neutron star cores. 
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7. Appendix: Numerical technics 

The integrals presented in this article, which were to be 
evaluated numerically, are of quite standard type in the field 
of degenerate Fermi fluids. However, as they involve various 
mathematical features that can easily lead to numerical 
difficulties, we have estimated that it was worth interacting 
with some applied mathematicians expert in sophisticated 
algorithms to dynamically control the accuracy of numerical 
multi-dimensional integrals. These people, J.-M. Chesneaux, 
F. Jezequel, F. Rico and M. Charrikhi from LIP6 of Paris 
6 University, are part of the "Cadna team", Cadna being 
the acronym for "Control of Accuracy and Debugging for 
Numerical Applications" and also the name of the numerical 
tool they created ('http://www-anp.lip6.fr/cadna/'). We 
shall not give here a description of their algorithm, but let us 
just mention that it is based on discrete stochastic arithmetic, 
and allows for calculations with a chosen precision (see, e.g., 
Chesneaux & Jezequel 119981 and Jezequel [2004 1 . 

Yet, the drawback of this method is that the time needed 
to carry out the computation can be very long, which prohibits 
to directly use it in physical simulations. As a consequence, we 
will no longer mention the Cadna code in the following, but 
it shall be noticed that we used it to test and calibrate our 
less sophisticated algorithms, and then to get clear ideas on 
the number of numerical points needed for various values of 
the parameters (see Jezequel et al. 2005 and Charrikhi et al. 
2002). Indeed, as we shall see now, the features of the integrals 
make that the needed number of points is not uniform in the 
(v = A/fc B T,£ = 8^i/k B T) plane. 

The main difficulties in the integrations to perform are 

- infinite size of the domains of integration; 

- singularities in the denominators of some factors (e.g. e x — 1 

for x — 0); 

- external free parameters (dimensionless gap amplitude v and 

dimensionless departure from equilibrium £) that can be 
very large or different. 

For not too precise calculations, the first of these issues 
is in fact quite easily dealt with, as using the asymptotic 
behaviours of exponential functions, one can show that it is 
not worth doing the integration on the full domain. Hence, the 
numerical domains are cut to some maximal values, which were 
chosen in such a way that exponential functions make further 
integration irrelevant. More precisely, due to the Fermi-Dirac 
distributions and to the "G-function" defined in Eq. (l22L one 
can see that the typical relevant scale for the "i-variable" of a 
superfluid nucleon is x ~ s x = £ + \/2v and that the scale 
for the neutrino variable is x v ~ s„ = £ + v. The maximal 
values for the numerical integration were 10 times the sum 
of 10 and of those typical scales, to prevent the occurence of 
problems for vanishing v and £. 

As far as singularities were concerned, we just dealt with 
them in putting in the code not exactly the functions but, 
beyond threshold values, some asymptotic formulae that are 
very good due to the exponential behaviours. The only trouble 



in this is that the presence of the possibly very large values 
of the v and £ parameters can make the cut-off values very 
huge. However, the total number of numerical points needed 
was kept quite reasonable using 

- some "logarithmic variables" (a suggestion done by V. 

Bezchastnov to K. Levenfish who kindly shared it with us): 
instead of directly integrating over the ^-variables, we made 
the integrations over some t-variables defined as 

x = s x (e* - 1) , (32) 

where s x is the typical scale defined above; 

- spectral decomposition: we evaluated the integrals with the 

so-called Gauss-Legendre quadrature method (see, e.g., 
Krvlov rrap^l . 

In this way, the main reason for needing huge num- 
bers of numerical points were either an anisotropic gap 
with a very large amplitude, or small gap amplitude 
but combined with large departure from beta equilib- 
rium. Anyhow, the subroutines used to produce the 
results displayed in Sect. |^| are gathered on the website 
1 http://luth2.obspm.fr/~etu/villain/Micro/Reduction.html ' , 
together with some comments and tests, and they include the 
determination of reasonable numbers of numerical points for a 
wide range of the v and £ parameters. Typically the proposed 
values (that can easily be changed playing in the subroutines) 
allow for a relative precision always at least better than 10 -3 
for v in [0, 10 3 ] and £ in [0, 10 4 ] . Of course, the subroutines 
can reach a much better precision just with a few additional 
numerical points and without huge change in the CPU time, 
but since there was no need for it, we limited the precision at 
that level. 

References 

Amundsen L., 0stgaard E., 1985, Nucl. Phys., A442, 163 
Andersson A., Comer G.L., 2001, Mon. Not. of the Royal 

Astron. Soc, 328, 1129 
Charrikhi M., Chesneaux J.-M., Jezequel F., Rico F., 

Villain L., 2002, proceeding of the SCAN2002 conference, 

Paris (France), 23-27 September 2002 
Chesneaux J.-M., Jezequel F., 1998, J. of Univ. Comp. Science, 

4(1), 2 

Fernandez R., Reisenegger, A., 2005, ApJ, 625, 291 
Friman B.L., Maxwell O.V., 1979, ApJ, 232, 541 
Gourgoulhon E., Haensel P., 1993, A&A, 271, 187 
Gusakov M.E., Haensel P., 2005, Nucl. Phys. A, 761, 333 
Haensel P., 1992, A&A, 262, 131 

Haensel P., Levenfish K.P, Yakovlev D.G., 2000, A&A, 357, 
1157 

Haensel P., Levenfish K.P, Yakovlev D.G., 2001, A&A, 372, 
130 

Haensel P., Levenfish K.P, Yakovlev D.G., 2002, A&A, 381, 
1080 

Jezequel F., 2004, App. Num. Maths., 50(2), 147 

Jezequel F., Rico F., Chesneaux J.-M., Charrikhi M., 2005, 

Math, and Comp. in Simulation, to appear. 
Krylov, V. I., 1962, Approximate calculation of integrals, ACM 

Monograph Series, The MacMillan Company, New York, 

NY 

Lattimer J.M., Pethick C.J., Prakash M., Haensel P., 1991, 
Phys. Rev. Lett., 66, 2701 



L. Villain & P. Haensel: Non-equilibrium beta processes in superfluid NSs 



11 



Lee U., 1995, A&A, 303, 515 

Levenfish K.P., Yakovlev D.G., 1994, Astron. Lett., 20, 43 
Lindblom L., Mendel G., 1995, ApJ, 444, 804 
Lindblom L., Owen B.J., Morsink S., 1998, Phys. Rev. Lett., 
80, 4843 

Lombardo U., Schulze H.-J., 2001, in: Physics of Neutron Star 
Interiors, D. Blaschke, N.K. Glendenning, A. Sedrakian 
(Eds.) Springer, Heidelberg, p. 30 

Prix R., Rieutord M., 2002, A&A, 393, 949 

Reisenegger A., Goldreich P., ApJ, 395, 240 

Reisenegger, A., 1995, ApJ, 442, 749 

Reisenegger, A., 1997, ApJ, 485, 313 

Shapiro S.L., & Teukolsky S.A., 1983, Black Holes, White 
Dwarfs and Neutron Stars, (Wiley-Interscience, New- York) 
Urpin V.A., Shalybkov D.A., 1996, MNRAS, 281, 145 
Villain L., Bonazzola S., Haensel P., 2005, Phys. Rev. D, 71, 
083001 

Yakovlev D.G., Levenfish K.P., 1995, A&A, 297, 717 
Yakovlev D.G., Levenfish K.P., Shibanov Yu.A., 1999, Physics- 

Uspekhi, 169, 825 
Yakovlev, D. G., Kaminker, A. D., Gnedin, O. Y., & Haensel, 

P., 2001, Phys. Rep. 354, 1. 



